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ABSTRACT 

The caustic technique uses galaxy redshifts alone to measure the escape velocity and 



^ ' mass profiles of galaxy clusters to clustrocentric distances well beyond the virial ra- 

dius, where dynamical equilibrium does not necessarily hold. We provide a detailed 
description of this technique and analyse its possible systematic errors. We apply the 
caustic technique to clusters with mass M200 f0 14 /i _1 Mq extracted from a cosmo- 
logical hydrodynamic simulation of a ACDM universe. With a few tens of redshifts 
per squared comoving megaparsec within the cluster, the caustic technique, on aver- 
age, recovers the profile of the escape velocity from the cluster with better than 10 
. percent accuracy up to r ~ 4r2oo- The caustic technique also recovers the mass profile 

with better than 10 percent accuracy in the range (0.6 — 4) r2oo> but it overestimates 
the mass up to 70 percent at smaller radii. This overestimate is a consequence of ne- 
glecting the radial dependence of the filling function J^(r). The l-er uncertainty on 
individual escape velocity profiles increases from ~ 20 to ~ 50 percent when the radius 
increases from r ~ 0.1r2oo to ~ 4r2oo- Individual mass profiles have 1-a uncertainty 
between 40 and 80 percent within the radial range (0.6 — 4)r2oo- When the correct 
KJi \ virial mass is known, the 1-a uncertainty reduces to a constant 50 percent on the 

same radial range. We show that the amplitude of these uncertainties is completely 
' due to the assumption of spherical symmetry, which is difficult to drop. Other poten- 

tial refinements of the technique are not crucial. We conclude that, when applied to 
individual clusters, the caustic technique generally provides accurate escape velocity 
and mass profiles, although, in some cases, the deviation from the real profile can be 
substantial. Alternatively, we can apply the technique to synthetic clusters obtained 
by stacking individual clusters: in this case, the 1-a uncertainty on the escape velocity 
profile is smaller than 20 percent out to 4r2oo- The caustic technique thus provides 
reliable average profiles which extend to regions difficult or impossible to probe with 
other techniques. 

Key words: gravitation - galaxies: clusters: general - techniques: miscellaneous - 
cosmology: miscellaneous - cosmology: dark matter - cosmology: large-scale structure 
of Universe 
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Clusters of galaxies are valuable tools to measure the cosmo- 
logical parameters and test str ucture formation models (e.g. 
IVoit| [2005: Diaf erio et al.ll2008l) and the galaxy-environment 
connection (e.g. Skibba et al.feodi^lHuertas-Companv et al.l 
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120091 ; iMartfnez etHI 120081 ; iDommguez et all l200ll ). The 
evolution of the cluster abundance is a sensitive probe of 
the cosmological parameters because clusters populate the 
exponential tail of the mass function of virialised galaxy sys- 
tems. Accurate mass measurements are however required to 
avoid the propagation of systematic errors into the estima- 
tion of the cosmological parameters. There are two families 
of mass estimators: those which estimate the mass profiles 
and those that measure the mass enclosed within a specific 
projected radius. 

Traditionally, the estimation of the cluster mass is based 
on the assumptions of spherical symmetry and dynamical 
equilibrium: either t he cluster gal axies move accordingly to 
the virial theorem (Z wickvl 119371 1 . or the hot intracluster 
plasma which emits in the X-ray band is in hydrostatic equi- 
libri um within th e gravitational potential well of the clus- 
ter (|Sarazinlll988l ). More sophisticated approaches apply the 
Jeans equation for a steady-state system, wit h the velocity 
anisotropy parame ter / 3 as a further u nknown |The fc White! 
ll98d : lMerrit"t]|l987l : see lBiviandliooll for a comprehensive re- 
view of various improvements of this method) . 

Dynamical equilibrium is also assumed when using the 
mass- X-ray temperature relation to esti mate the cluste r 
mass (e.g. iPierpaoli et~aH 120031 : see also iBorganil [2006). 
Both the slope and the normalization of the observed mass- 
temperature relation is grossly reproduced by synthetic 
clusters obtained wit h JV-body/hydrodynamical simulations 
l|Borgani et al.ll2004 ). However, the presence of non-thermal 
pressure support (e.g., related to turbulent gas motions) and 
the complex thermal structure of the intra-cluster medium 
can significantly bias cluster mass estimates based on the 
application of hydrosta tic equilibrium to the X-ray esti- 
mated temperature (e.g . Rasia et a l. 2006; Nagai et al]|2007l ; 



IPiffaretti fc ValdarniniH2008h . In principle, one can improve 
the mass estimate by exploiting the integrated Sunyaev- 
Zel'dovich effect, which depends on the first power of the gas 
density, rather than the square of the density as in the X-ray 
emission, and yields a correlation with mass whi ch is tighter 
than the mass-X-r ay temperature relation (e.g. iMotl et al.l 
120051 : [Nagai 200(3). Only recently, such results have been 
confirmed by observat ions of real clusters l|Rines et al.ll201ol ; 
lAndersson et al]|2010h . However, their confirmation over a 
large statistical ensemble has to await the results from ongo- 
ing large Sunyaev-Zel'dovich surveys. Alternatively, one can 
bypass the problematic gas physics and use the correlation 
between mass and optical richness, which is relatively easy 
to obtain from observations in the optic al/near-IR bands, 
but r eturns mass with a poorer accuracy (jAndreon fc Hurnl 
l2010l ). 

The dynamical equilibrium assumption can be dropped 
in the gravitational lensing techniques, because the lensing 
effect depends only on the amount of mass along the line o f 
sight and not on its dynamical state (e.g. ISchne idcr 2006). 
It is relevant to emphasize that all these methods do not 
measure the cluster mass on the same scale: optical obser- 
vations measure the mass within ~ V2oo, where r2oo is the 
radius within which the average mass density is 200 times 
the critical density of the universe; X-ray estimates rarely 
go beyond ~ 0.5r2oo, where the X-ray surface brightness be- 
comes smaller than the X-ray telescope sensitivity; lensing 
measures the central mass within ~ 0.1r2oo or the outer re- 
gions at radii larger than ~ V2oo depending on whether the 



strong or weak regime applies. Scaling relations do not pro- 
vide any information on the mass profile and give the total 
mass within a radius depending on the scaling relation used, 
but still within ~ T2oo- 

iDiaferio fc Gelled (|l997l ) (DG97, hereafter) suggested a 
novel method, the caustic technique, to estimate the mass 
from the central region out to well beyond r2oo with galaxy 
celestial coordinates and reds hifts alone an d without as- 
suming dynamical equilibr ium (lDiaferidl2009|). Prompted by 
the iV -body simulations of lvan Haarlem fc van de Weygaert] 
(1993), DG97 noticed that in hierarchical models of struc- 
ture formation, the velocity field in the cluster outskirts is 
not perfectly radial, as expected in th e spherical infall model 
l|Regos fc Gelled Il989l : iHiotelid l200ll ) but has a substantial 
random component. They thus suggested to exploit this fact 
to extract the galaxy escape velocities as a function of ra- 
dius from the distribution of galaxies in redshift space. In the 
caustic method, the velocity anisotropy parameter /3 and the 
mass density profile in hierarchical clustering models com- 
bine in such a way that their knowledge is largely unneces- 
sary in estimating the mass profile. This property explains 
the power of the method. 

This method is particularly relevant because it is an 
alternative to lensing to measure mass in the cluster external 
regions and, unlike lensing, it can be applied to clusters at 
any redshift, provided there are enough galaxies to sample 
the redshift diagram properly. We will see below that a few 
tens of redshifts per squared comoving megaparsec within 
the cluster are sufficient to apply the technique reliably. This 
request might have appeared demanding a decade ago, but 
it is perfectly feasible for the large redshift surveys currently 
available. 

iGeller. Diaferio. fc Kurtzl (QUI) were the first to ap- 
ply the caustic method: they measured the mass profile 
of Coma out to 10ft _1 Mpc fro m the cluster centre and 
were a ble to demonstrate that the lNavarro. Frenk. fc White] 
l|l997h (NFW) profile fits well the cluster density profile 
out to these very large radii, thus ruling out the isother- 
mal sphere as a viable cluster model; a few years later, the 
failure of the isothermal model was confirmed by the first 
analysis based on gravitational lensing applied to CI 0024 
j|Kneib et all 120031 '). The goodness of the NFW fit out to 
5- Wh- 1 Mpc was confirmed by applying the caustic tech- 
nique to a sample of nine clusters densely sampled in their 
outer regi ons, the Cluster A nd Infall Region Nearby Survey 
(CAIRNS, Rincs et al. 200j|), and later to a complete sample 
of 72 X-ray selected clusters with galaxy redshifts extracted 
from the Fourth Data Release of the Sloan Digital Sky Sur- 
vey (C luster Infall Regions in the Sloan Digital Sky Survey: 



CIRS, Rincs & Diaferio 



2006) and from the Fifth Data Re- 



lease (|Rines fc Diaferid l2010). CIRS is currently the largest 



sample of cluste rs whose mass pr o files h ave been measured 
out to ~ 3r2oo; iRines fc Diaferio] (I2006T ) were thus able to 
obtain a statistically significant estimate of the ratio be- 
tween the masses within the infall and the virial regions: 
they found a value of 2.2 ± 0.2, in agreement with current 
models of cluster formation. These analy ses have been ex- 
tende d to a sample of groups of galaxies (|Rines fc Diaferid 
l20ld) . IRines. Diaferio. fc Nataraianl (l2007n ~ also used the 
CIRS sample to estimate the virial mass function of nearby 
clusters and determi ned cosmological par ameters consistent 
with WMAP values l|Dunklev et al.ll2009h . 
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A good fit with the NFW profile o ut to ~ 2r2oo was 
also found by iBiviano fc Girardil <|2003h who applied the 
caustic technique to an ensemble cluster obtained by stack- 
ing 43 cluster s from the Two De gree Galaxy Redshift Sur- 
vey (2dGFRS, IColless et aLlbOOll ). Here, unlike the previous 
analyses, the caustic method was not applied to individual 
clusters because the n umber of galaxies p er cluster was rela- 
tively small. Recently. iLemze et al.l (|2009l ) have applied both 
the caustic technique and the Jeans analysis to ~ 500 clus- 
ter members of A1689. The estimated virial mass from both 
methods agrees with previous lensing and X-ray analyses. 

The caustic method does not rely on the dy- 
namical state of the cluster and of its external re- 
gions: there are therefore estimates of the mass of 
unrelaxed systems, namely the Shapley supercluster 
l|Reisenegger et al ] |2000| ; iProust et aT] l2006to . the Fornax 
poor cluster, which shows two distinct d ynamical compo- 
nents (iDrinkwater. Gregg, fc Collessll200ll ). the A2199 com- 
plex (|Rines et alj|2002h . and two cluster s, A168 and A1367 , 
which are und ergoing major merger s (|Rines et al.l [2003). 
More recently, iDiaferio et all (l2005bh analysed CI 0024, a 
cluster that is likely to have suffered a recent merging event 
|Czoske et al.1120021 ). and showed that the caustic mass and 
the lensing mass agree with each other, but disagree with the 
X-ray mass, which is the only estimate relying on dynamical 
equilibrium. 

The method was shown to return reliable mass profiles 
when applied to synthetic clusters extracted from iV-body 
simulations. However, these analyses of the performance of 
the method either were done when the operative procedure 
of the technique was not yet completely developed (DG97) 
or the results of a fully detailed analy sis of the systematic er- 
rors was not provided (|Diaferidll999T ) (D99, hereafter). The 
dense redshift surveys currently available, especially around 
clusters, provide ideal data sets where the caustic technique 
can now be applied. With this perspective, it is therefore 
timely to reconsider the possible systematic errors and bi- 
ases of this technique, in view of a robust application to the 
new data sets. 

The purpose of this paper is to provide both a transpar- 
ent description of the technique, with some improvements to 
what is described in D99, and a thorough statistical analy- 
sis. The caustic technique represents a powerful method to 
infer cluster mass profiles and complements other methods 
based on X-ray, Sunyaev-Zeldovich and lensing observations 
which probe different scales of the clusters. 

In Section [5] we describe the basic idea of the caustic 
technique. Section [3] describes the simulated cluster sam- 
ple used, and Section U describes the implementation of the 
technique. In Sections 15. 21 and 15. 31 we estimate the accuracy 
of the escape velocity and mass profiles of galaxy clusters 
estimated with the caustic technique. In Section[S]we inves- 
tigate the systematics due to the choice of the parameters. 
We summarize our main results and draw our conclusions 
in Section [7] 



2 BASICS 

In this section, we briefly review the simple physical idea 
behind the caustic technique. More details are given in DG97 
and D99. 



In hierarchical clustering models of structure formation, 
clusters form by the aggregation of smaller systems falling 
onto the cluster from the surrounding re gion. The accretio n 
does not take place purely radially (e.g. White et al.ll201oT ). 
therefore galaxies or dark matter particles within the falling 
clumps have velocities with a substantial non-radial compo- 
nent. Specifically, the r.m.s. {v 2 ) of these velocities is due 
to the gravitational potential of the cluster and the groups 
where the galaxy resides, and to the tidal fields of the sur- 
rounding region. When viewed in the redshift diagram, viz. 
the plane of the line-of-sight (l.o.s.) velocity w.r.t. the clus- 
ter centre and the clustrocentric radius r, galaxies populate 
a region with a characteristic trumpet shape with decreasing 
amplitude A with increasing r. This amplitude is related to 
{v 2 }. The breakthrough of DG97 was to identify this ampli- 
tude with the escape velocity from the cluster region cor- 
rected for a function depending on the velocity anisotropy 
parameter j3. 

Assuming a spherically symmetric system, the escape 
velocity v 2 sc (r) — —2cf>(r), where cj)( r ) ls the gravitational 
potential originated by the cluster, is a non-increasing func- 
tion of r, because gravity is always attractive and d(/>/dr > 0. 
At any given radius r, we expect that observing a galaxy 
with a velocity larger than the escape velocity is unlikely. 
Thus, we identify the escape velocity with the maximum ve- 
locity that can be observed. It follows that the amplitude 
A at the projected radius r± measures the average com- 
ponent along the l.o.s. of the escape velocity at the three- 
dimensional radius r = r±. To determine this average com- 
ponent of the velocity, we consider the velocity anisotropy 
parameter /3(r) = 1 — ((v$} + {v 2 ))/2{v 2 } , where vg, v^, and 
v r are the longitudinal, azimuthal and radial components 
of the velocity v of a galaxy, respectively, and the brack- 
ets indicate an average over the velocities of the galaxies in 
the volume d 3 r centred on position r. If the cluster rota- 
tion is negligible, {v 2 .} — {v 2 } — (vf os ), where vi os is the l.o.s. 
component of the velocity, and {v 2 } — {v 2 } — 2{v? OB ). By sub- 
stituting in the definition of /3, we obtain {v 2 } = (v 2 os )g((3) 
where 



2/3(r) 



/3(r) 



(1) 



By applying this relation to the escape velocity at radius r, 
{wicM) = -20(r), and assuming that A 2 {r) = (l4c,los}» we 
obtain the fundamental relation between the gravitational 
potential (j)(r) and the observable caustic amplitude A(r) 



2^{r)=A 2 {r)g{P) = <j ) p{r)g{P) 



(2) 



Note that the gravitational potential profile is related to the 
caustic amplitude by the function g(/3). Therefore, after the 
caustic amplitude estimation, /3 becomes the only unknown 
function for the gravitational potential estimation. 

The further novel suggestion of DG97 is to use this rela- 
tion to infer the cluster mass to very large radii. To do so, one 
first notices that the mass of an infinitesimal shell can be cast 
in the form Gdm = -20(r)J r (r) dr = A 2 {r)g(P)F{r) dr 
where 



T{r) = -2-kG 



p(r)r 2 



(3) 



Therefore the mass profile is 
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GM(< r) 



A (r)Fp(r) dr 



(4) 



where ^(r) = F{r)g{P). 

Equation Q however only relates the mass profile to 
the density profile of a spherical system and a profile can 
not be inferred without knowing the other. DG97 solve this 
impasse by noticing that in hierarchical clustering scenarios, 
J-(r) is not a strong function of r. This is easily seen in the 
case of the NFW model which is an excellent description 
of dark matter density profiles in the hierarchical clustering 



scenario: 



.^nfwM 



2(r + r s ) 2 ln(l + r/r s ) ' 



(5) 



where r s is a scale factor. One expects that Tp{r) is also 
a slowly changing function of r if g(/3) is. DG97 and D99 
show that this is the case, and here we confirm their result. 
The final, somewhat strong, assumption is then to consider 
J-p(r) = Tf) = const altogether and adopt the recipe 



GM{<r) =F P / A 2 (r)dr 



(6) 



It is appropriate to emphasize that equations ([2]) and 
@ are rigorously correct, whereas equation is a heuris- 
tic recipe to estimate the mass profile. Based on equation 
((2| the caustic technique can also estimate a combination 
of the gravitational potential profile and velocity anisotropy 
parameter /3. Below, we show the accuracy of these esti- 
mates. 



3 THE SIMULATED CLUSTER SAMPLE 

Before analysing the systematic errors of the caustic tech- 
nique, we need to check how well the physical assumptions 
of the method are satisfied. We do so by using TV-body sim- 
ulations, assuming that these simulations are a faithful rep- 
resentation of the real universe. 

Here we use the cluster sample extracted from the cos- 
mological iV-body/h ydrodynamical simulation described in 
iBorgani et all \20oS . They simulate a cubic volume, 192 
h- 1 Mpc on a side, of a flat ACDM universe, with mat- 
ter density flo = 0.3, Hubble parameter Ho = 100/t km 
s _1 Mpc -1 with h = 0.7, power spectrum normalization 
as — 0.8, and baryon density fib = 0.02/i -2 . The den- 
sity field is sampled with 480 3 dark matter particles and 
an initially equal number of gas particles, with masses 



771DM 



4.6 x 10 9 /i _1 M ra and m„ 



6.9 x 10 8 /7 _1 M, 



©, re- 



spectively. The Plummer-equivalent gravitational softening 
is set to 7.5 ft -1 physical kpc at z < 2, while being fixed in 
comoving units at higher redshift. 

The simulation wa s run with GADGET-2 

|Springel. Yoshida. fc White! 120011 ). a massively parallel 
Tree+SPH code with fully adaptive time-stepping. Besides 
standard hydrodynamics, it also contains a treatment of 
radiative gas cooling, star formation in a multi-phase inter- 
stellar medium and fee dback from supernovae in the form 
of gal actic ejecta (e.g., IBorgani et al.l 120041 ; iDiaferio et al.l 
l2005ah . Here we are only interested in the gravitational 
dynamics of the matter distribution. 

The simulation volume yields a cluster sample large 
enough for statistical purposes. We identify clusters in the 
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Figure 1. Concentration parameters c vs. the cluster mass 
M|qq W of our simulated cluster sample. 




Figure 2. Profiles of the ratio between the mass profile of each 
cluster predicted by the NFW fit and its true mass profile: 90 
(50) percent of the profiles are within the upper and lower dotted 
(solid) curves. The central solid curve is the median profile. For 
each cluster, the NFW fit is only performed to the mass distribu- 
tion within l/i" 1 Mpc. 



simulation box with a two-step procedure: a friends-of- 
friends algorithm applied to the dark matter particles alone 
provides a list of halos whose centres are used as input to 
the spherical ov erdensity algorithm which outputs the final 
list of clusters (jBorgani et al.l 120041 ). Centred on the most 
bound particle of each cluster, the sphere with virial over- 
density 200, with respect to the critical density, defines the 
virial radius 7-200- At redshift z — the simulation box con- 
tains 100 clusters with mass M(< r 2 oo) > 1O 14 /i _1 M . This 
cluster set is our reference sample. 



3.1 Cluster properties in 3D 

We fit the three-dimensional (3D) real cumulative mass pro- 
file with the NFW model 



M(< r) = M s 



In 1 + - - 



r/r s 



1 + r/r s 



(7) 
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Figure 3. Profiles of the ratio between the gravitational potential 
profile of each cluster predicted by the NFW fit and the numerical 
gravitational potential profile derived from the true mass distri- 
bution within r max = 10h~ 1 Mpc from the cluster centre: 90 
(50) percent of the profiles are within the upper and lower dot- 
ted (solid) curves. The central solid curve is the median profile. 
The lower and upper dashed curves are the median profiles for 
r ma i = 8 and 6/i _1 Mpc, respectively. 



where M s = M(< r s )/(ln2 - 1/2) = AnS c p CI r 3 s . 
(200/3)c 3 /[ln(l + c) - c/(l + c)], and c = r 2 N F W /r s . 



5c = 
We fit 



the mass profile rather than the density prof ile to resemble 
the pr ocedure adopted with real clusters (e.g. iDiaferio et al.l 
l2005bl ). Similarly, for the fit, we only consider the mass dis- 
tribution within r = l/i -1 Mpc. Moreover, the parameters 
Ms and r s are less correlated than 8 C and r a and adopting 
them in the fit provides more robust results (|Mahdavi et al.l 
I1999T ). 

From the NFW fit we derive the parameters 
r M0 W an d c The relation between c and Af|Jo W = 
(47r/3)(r^ t 1 o W ) 3 200p cr is shown in Figure Q] Our clus- 
ter sample provides the percentile ranges rgoo = 
[0.77, 0.86, 1.17]7i _1 Mpc and c = [2.06, 4.62, 6.50]Q The ra- 
dius r2oo W is basically identical to the true r2oo derived from 
the actual mass distribution. Their ratio is r^)o W /r2oo = 
[0.99, 1.00, 1.02]. The NFW model is an excellent fit to the 
true mass profile (see Figure [2]) within rzw- Therefore the 
ratio between the mass M^txt™ and the actual M200 is corre- 
spondingly close to 1: M 2 ^ w /M 20 o = [0.96, 1.00, 1.05]. The 
NFW fit is on average very good also at radii larger than 
r 200, although the scatter substantially increases. In the very 
centre, the NFW model underestimates the mass by 50 per- 
cent. This is due to the fact that the total mass includes dark 
matter, gas, and stars and in the simulation there is always 
a central galaxy which is unrealistically massive: the stellar 
particles, on average, contribute 50 (15) percent of the total 
mass within 0.02 (0.1) r2oo- When the NFW fitting routine 
tries to accomodate the mass profile within l/i -1 Mpc, un- 
derestimating the mass profile in this very central region 
yields the minimum \ 2 • 



1 Throughout this paper the notation [91,92,93] shows the me- 
dian 92 and the range [91,93] which contains 90 percent of the 
sample. 



The caustic amplitude is related to the gravitational 
potential <j>(r) and we are thus interested in determining 
(j)(r) in our simulation. For an isolated spherical system with 
density profile p(r), the potential ci(r) obeying the Poisson 
equation can be cast in the form 



(r) 



-4ttG 



p(x)x dx + p(x)xdx 



(8) 



In the real universe, the relevant quantity is the gravita- 
tional potential originating from the mass density fluctua- 
tions around the mean density (p). We can thus use equa- 
tion (|8]) for non-isolated clusters by replacing p(r) with 
p(r) — (p) = (p)S(r). The second integral is finite, because at 
sufficiently large r, S(r) ~ 0. In the simulation, we replace 
the upper limit of the second integral with r max = lO/i -1 
Mpc. Figure [3] compares (!> nnm (r) computed from the actual 
mass distribution around each cluster in the simulation, with 
the gravitational potential 



<j>NFw(r) = 



r \ r 3 



(9) 



expected from an isolated cluster described by the NFW 
model. The NFW potential well is deeper than the numer- 
ical 4>num(r), because it neglects the mass surrounding the 
cluster. This mass exerts a pull to the mass within the clus- 
ter that makes the actual <^ n um(r) to be 10-30 percent shal- 
lower. The upper limit r max of the second integral of equa- 
tion JHJ, adopted in the numerical estimate, plays a neg- 
ligible effect when it is large enough: the median <f) n um(r) 
for rmax = 8/t -1 Mpc is indistinguishable from the profile 
computed with r max = 10/i _1 Mpc. 

Figure [3] shows the profiles of the velocity anisotropy 
parameter j3(r) and the other functions defined in Section [2] 
for our sample of 100 simulated galaxy clusters. This figure 
confirms the results of the ACD M model presented in Figure 
3 of D99 (see also Figure 25 in IDiaferio et al.ll200ll . which 
shows the velocity field of simulated galaxies rather than 
dark matter particles). Specifically, Figure [3] shows that, on 
average, /3(r) < 0.7 at r < 3r2oo and consequently, on these 
scales, the median profile of g((3) varies by less than 30 per- 
cent. Individual clusters may of course have wider variations. 

We remind that fi(r) is derived from the velocities of 
all the particles in the simulations (dark matter, gas and 
stars). These velocities are set in the rest frame of each 
cluster and are augmented of the Hubble flow contribution 
Hor. This term provides an increasing contribution to the 
particle velocity which is not negligible at r > (1 — 2)r2oo- 
Specifically, by assuming dynamical equilibrium and setting 
v pi r ) ~ GM(< r)/r, where v p (r) is the average peculiar 
velocity at r, we see that v p decreases with r; therefore, 
the total velocity v(r) = v p (r-) + Hor reaches a minimum 
value when v p (r) ~ H Q r, viz. when r ~ [GM (< r) /Ho] 1/:i ~ 
4.01ft" 1 Mpc, where we have set M = 1.49 x 10 14 h~ 1 M Q , 
the median Af 2 oo of our cluster sample. Since the tangential 
component of the total velocity is unaffected by the Hubble 
flow, the minimum of v p (r), and consequently of (fj:), is a 
minimum of f3. For our sample the median r 2 oo = 0.86/t" 1 
Mpc, and we have a minimum /? at r/r 2 oo = 4.66 which is 
in rough agreement with Figure U 

Finally, Figure U shows that J-(r) and J-p{r) are slowly 
varying functions of r, as expected. We discuss these profiles 
in Section \5. 31 
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Figure 4. Profiles of the functions /3(r), g(0), T(r) and Ta(r) described in the text: 50, 68 and 90 percent of the profiles are within the 
upper and lower dotted, solid and dashed curves. The central solid curves are the median profiles. 



3.2 The mock cluster catalogue 

We compile mock redshift catalogues from our 100 simu- 
lated clusters as follows. We project each cluster along ten 
random lines of sight; for each line of sight, we also project 
the clusters along two additional lines of sight in order to 
have a set of three orthogonal projections for each random 
direction. We thus have a total of 100 x 10 x 3 = 3000 mock 
redshift surveys. We locate the cluster centre at the celes- 
tial coordinate (a, 8) = (6 h ,0°) and redshift cz = 32000 
km s . We consider a field of view of 12ft" 1 Mpc on a side 
at the cluster distance. The simulation box is 192ft _1 Mpc 
on a side, and the mock survey contains a large fraction of 
foreground and background large-scale structure. We only 
consider a random sample of 1000 dark matter particles in 
each mock catalogue. Only a fraction of these 1000 parti- 



cles are within 3r2oo of the cluster centre, specifically this 
number has percentile range [96, 185, 408] in our 3000 mock 
catalogues. In Sect. l6~4l we investigate the effect of varying 
the number of particles in the mock survey. 



4 THE CAUSTIC TECHNIQUE 

The celestial coordinates and redshifts of the galaxies within 
the cluster field of view are the input data of the caustic 
technique. The technique proceeds along four major steps: 

(1) arrange the galaxies in a binary tree according to a 
hierarchical method; 

(2) select two thresholds to cut the tree twice, at an up- 
per and a lower level: the largest group obtained from the 
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upper-level threshold identifies the cluster candidate mem- 
bers; the lower-level threshold identifies the substructures 
of the cluster; the cluster candidate members determine the 
cluster centre, its velocity dispersion and size; 

(3) with the cluster centre of the candidate members, 
build the redshift diagram of all the galaxies in the field; with 
the velocity dispersion and size of the candidate members 
determine the threshold k which enters the caustic equation, 
locate the caustics, and identify the final cluster members; 

(4) the caustic amplitude determines the escape velocity 
and mass profiles. 

We detail these steps in turn. 



4.1 Binary tree 

Each galaxy is located by the vector n — (on,6i,ri), where 
Qi and 5i are the celestial coordinates, and rt is the comov- 
ing distance to the observer of a source at redshift Zi, i.e. 
the spatial part of the geodesic travelled by a light signal 
from the source to the observer. The comoving radial coor- 
dinate rj is determined by the relation J r * dx/y/1 — kx 2 = 
J Q Zl cdz/H(z), where k — (Ho / c) 2 (fio + fiA — 1) is the space 
curvature and H 2 (z) = H 2 [fi (l+2) 3 +(l-fio-fiA)(l+.z) 2 + 
Qa] is the Hubble constant at redshift z. In the Einstein-de 
Sitter universe (k = 0, fio — 1) 



_2c 

Ho 



(1 + %) 1/2 



(10) 



The comoving separation ria between two sources with co- 
moving radial coordinates r\ and r 2 and angular separation 
9 as observed from O is 



2 
1-12 



2 , 2 2 2/1, 2 

r-± + r 2 — «r 1 r 2 (l + cos 



2r\r 2 cos9x 



1 — ttrf 



1 — «r§ 



(11) 



This equation is a generalization of th e cosine law o n the 
two-dimensional surface of a sphere (e. g. Peebles 1993, equa- 
tion 12.42). 

For the construction of the binary tree, we are interested 
in defining an appropriate similarity based on the pairwise 
separation m. In general, the rank of the separation ri2 of a 
number of sources with different redshift z% is not preserved 
when we vary fio and fiA- However, a sufficient condition 
for having the ranks preserved is that the universe is never 
collapsing and k $S 0. This condition is satisfied when fio 1 
and fio + fiA ^ 1. 

Based on this argument, we can compute ria in the 
Einstein-de Sitter universe, and use the ordinary Euclidean 
geometry to derive the binary tree similarity. For each galaxy 
pair, we define the l.o.s. vector I = (n + ra)/2 and the 
separation vector ri2 = T\ — r%. We then determine the 
component of ri2 along the line of sight, 



ri2 



I 



\l\ 



2 2 

\ri +r 2 \ 



(12) 



and the component perpendicular to the line of sight 

r P = (r? 2 -^ 2 ) 1/2 . (13) 

We now consider the proper l.o.s. velocity difference n = 
7r/(l + Zi) and the proper projected spatial separation R p = 
r p /(l + zi) of the galaxy pair at the intermediate redshift 



zi that satisfies the relation n 
similarity 



(n + r 2 )/2. We adopt the 



E, 



^mirrii 1 rriimi 
= —G „ + -- 



2 rrii + rrij 



-a 2 . (14) 



This similarity reduces to the one adopted by D99 at small 
2. In equation (|14[) . the galaxy mass is a free parameter. 
D99 assumes m, = rrij = 10 12 /i _1 Mq, as we do here. 
Varying the galaxy mass changes the details of the binary 
tree, because one differently weights n and R p . However, in 
the range (10 11 , 10 13 )/i _1 Mq, the centre location and the 
galaxy membership remain substantially unchanged. One 
can of course include different galaxy masses proportional 
to the galaxy luminosities. This prescription turns out to be 
unnecessary, because it would not increase the accuracy of 
the determination of the mass and escape velocity profiles. 
In fact, we will see below (Section 16. 3[) that projection ef- 
fects completely dominate the uncertainties of the caustic 
method. 

For the sake of completeness, we note that, in a non- 
flat universe, with fio < 1 and fio + fiA < 1, one can still 
perform the component separation with equations (|12[) and 
(|13[1 on a flat space tangent to the observer's location point 
O, provided one defines an appropriate bi-unique function 
between a geodesic distance r* and a vector n on the flat 
space. 

The binary tree is built as follows: (i) each galaxy can 
be thought as a group g a \ (ii) to each group pair g a ,gp it is 
associated the binding energy E a p — min{_Eij}, where iJy 
is the binding energy between the galaxy i £ g a and the 
galaxy j £ gp; (iii) the two groups with the smallest binding 
energy E a p are replaced with a single group g 7 and the total 
number of groups is decreased by one; (iv) the procedure is 
repeated from step (ii) until we are left with only one group. 
Figure [5] shows the binary tree of a simulated cluster as an 
example. 



4.2 Cutting the binary tree: The a plateau 

Gravity is a long-range force and defining a system of galax- 
ies is somewhat arbitrary. The most widely used approach is 
to select a system according to a density threshold, either in 
redshift space, for real data, or in real space, for simulated 
data. 

When applied to real data, the binary tree described in 
the previous section has the advantage of building a hierar- 
chy of galaxy pairs with increasing, albeit projected, binding 
energy. The tree automatically arranges the galaxies in po- 
tentially distinct groups based on a single parameter, the 
galaxy mass. To get effectively distinct groups, however, we 
need to choose a threshold, i.e. the node from which the 
group members hang, where we cut the tree. This choice has 
been traditionally arbitrary l|Materne]| 19781 ; ISerna fc Gerball 
Il996h . 

D99 suggests a more objective criterion based on the fol- 
lowing argument. We first identify the main branch of the bi- 
nary tree, that is the branch that links the nodes from which, 
at each level of the tree, the largest number of leaves hang. 
The velocity dispersion a x of the galaxies hanging from a 
given node x shows a characteristic behaviour when walking 
towards the leaves along the main branch (Figure [6j : ini- 
tially it decreases rapidly, it then reaches a long plateau and 
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Figure 5. Dendrogram representation of the binary tree of a subsample of the particles in the field of a simulated cluster. The particles 
are the leaves of the tree at the bottom of the plot. The particles within 3r200 in real space are highlighted in black. The thick path 
highlights the main branch. The horizontal lines show the upper and lower thresholds used to cut the tree. Only as a guide, some nodes 
are labelled on the left-hand side, with their number of descendants in brackets. 



again drops rapidly towards the end of the walk. The plateau 
is a clear indication of the presence of the nearly isothermal 
cluster: at the beginning of the walk, a x is large because of 
the presence of background and foreground galaxies; at the 
end of the walk the tree splits into different substructures 
and a x drops again. 

The two nodes xi and X2 that limit the a plateau are 
good candidates for the cluster/substructure identification. 
To locate x\ and x%, we proceed as follows. We consider 
the density distribution of the a x . The mode of this density 
distribution is a p \. We then identify the Ns nodes whose 
velocity dispersion a x fulfills the inequality 



|o" p i — o x 

CTpl 



s£ 5 



(15) 



Clearly, the number of nodes Ns depends on the param- 
eter 5. We limit the parameter 8 in the range [0.03, 0.1], but 
we also compute the number of nodes N0.3 when 8 = 0.3. 
We determine the number of nodes Ns corresponding to in- 
creasing values of 5 £ [0.03,0.1] by step of 0.01 until Ns is 
larger than O.8-/V0.3. This sets the final number of nodes Ns. 



If all the N s with 8 £ [0.03,0.1] are smaller than 0.8iVo.3, 
we choose O.8-/V0.3 as the final number of nodes. The up- 
per limit of the range of 8, 8 = 0.1, is chosen because it 
always provides a sufficiently large number of nodes (15 in 
our sample); the arbitrary threshold O.8-/V0.3 is chosen be- 
cause it enables to deal efficiently with particularly peaked 
density distributions of a x . Among the Ns nodes, we locate 
the five nodes closest to the root and the five nodes closest 
to the leaves; among the former set, the final node x\ has a x 
with the smallest discrepancy from cr p \, and similarly for the 
final node X2 among the five nodes closest to the leaves. This 
procedure, illustrated in Figure [7] guarantees a sufficiently 
large number of nodes along the a plateau and prevents us 
from locating the extrema of the a plateau on nodes whose 
a x are too discrepant from <j p i. 
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Figure 7. Line-of-sight velocity dispersion of the leaves of each node along the main branch of the binary tree shown in Figure [5] The 
filled red circles are the Ng nodes. The dashed line indicates cr p i, whereas the dot-dashed lines show the range denned by equation H15I I: 
8 = 0.05 in this case. The ten circles with a yellow centre indicate the five nodes closest to the root and the five nodes closest to the 
leaves. These nodes also appear in the two insets. The two blue squares denote the two final nodes xi and X2- 




012345 1 2 3 4 

r (Mpc/h) r (Mpc/h) 



Figure 8. Left panel: the function f q (r,v) on the redshift diagram (colour map and contours) of the cluster shown in Figures I5I7I The 
thick line shows where f q (r,v) = k. Right panel: the redshift diagram of the same cluster with the upper and the lower caustics. The 
black and cyan lines are the estimated and true caustics respectively; the error bars are estimated with the procedure described in section 
15.51 The symbols are the particles in the catalogue: the red dots are the particles within a sphere of 3r200 centred on the cluster centre, 
in real space; the green diamonds are the particles whose l.o.s. velocity exceeds 3.5 times the velocity dispersion (t; 2 ) 1 / 2 of the candidate 
members. 
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Figure 6. Velocity dispersion of the leaves of each node along 
the main branch of the binary tree shown in Figure [5] The filled 
square and triangle indicate the final thresholds (nodes x\ and 
X2) chosen by the algorithm. The dashed line shows the l.o.s. 
velocity dispersion of the particles within a 3D sphere of radius 
3r-200- 



4.3 Redshift-space diagram, caustics and final 
members 

The candidate cluster members are the galaxies hanging 
from node xi and constitute the main group of the binary 
tree. These galaxies determine the centre and the size R 
of the cluster. The redshift z c of the cluster is the me- 
dian of the candidate redshifts. The celestial coordinates 
of the centre are the coordinates of the two-dimensional 
(2D) density peak of the candidates. To find the peak, we 
compute the 2D density distribution f q (a,5) of the can- 
didates on the sky with the kernel method described be- 
low (equation I18[) . In general, the smoothing parameter 
h c is automatically chosen by the adaptive kernel method. 
However, to save substantial computing time, here we set 
h c = 0.15(Da/320/i -1 Mpc) rad, where Da(z c ) is the an- 
gular diameter distance of the cluster. This choice yields 
accurate results. The cluster size R is the mean projected 
separation of the candidate cluster members from the cluster 
centre. 

We then build the redshift diagram, which is the plane 
(r, v) of the projected distance r and the l.o.s. velocity v of 
the galaxies from the cluster centre. If tp is the angular sep- 
aration between the cluster centre and a galaxy at redshift 
z, then 



cD A (z c ) . . 
r — — ^rr — - sin tp , 



and 



Ho 



l + z c 



(16) 



(17) 



where we have assumed that the galaxy velocity within the 
cluster is much smaller than the speed of light c and we have 
neglected the peculiar velocity of the cluster. Note that to 
avoid artificial depletion of the caustic amplitude at small r 
due to the small number of galaxies in the central region, the 
galaxy distribution is mirrored to negative r. As an exam- 



ple, in Figure [8] we show the redshift diagram of the system 
shown in Figures [5jT] 

It is now necessary to locate the caustics. The distribu- 
tion of N galaxies in the redshift diagram is described by 
the 2D distribution: 



-f-K 

N^h 2 



x, 



hi 



where x = (r, v) 



the adaptive kernel is 



K(t) 



t" 







t < 1 

otherwise; 



(18) 



(19) 



hi — h c h op tXi is a local smoothing length, Aj = 
[7//i(Xi)] 1/,2 ; fi is equation (fT8|) where h c = Xi = 1 for 
any i, and log 7 = £^ log[/i(xj)]/jV. D99 describes how the 
adaptive kernel method automatically determines h c . 
The optimal smoothing length 



hop 



6.24 

ATVe 



„2 , „2\ 1/2 



(20) 



depends on the marginal standard deviations o r and a v of 
the galaxy coordinates in the (r, v) plane. These two cr's 
must have the same units and we require the smoothing 
to mirror the uncertainty in the determination of the posi- 
tion and velocity of the galaxies. We therefore rescale the 
coordinates such that q = a v ja r assumes a fixed value. We 
choose q — 25. With this value, an uncertainty of 100 km s -1 
in the v direction, for example, weights like an uncertainty 
of 0.02ft -1 Mpc in the r direction. The optimal smoothing 
length /i op t adopted by D99 is half the value we use here. 
However, that h opt is derived in the context of the proba- 
bility density estimate under the assumption that x = (r, v) 
is a normally dis tributed random variate with unit variance 
l)Silvermanll 19861 . sects. 4.3 and 5.3). However, this assump- 
tion clearly does not apply to our context here, and we find 
that we obtain more accurate results with equation (|20p . 

The caustics are now the loci of the pairs (r,v) that 
satisfy the caustic equation 



(21) 



where k is a parameter we determine below. In general, the 
density distribution f q (r,v) is not symmetric around the 
straight line v — in the (r, v) plane, and equation (|21[) de- 
termines two curves v up (r) and Vdown(r), above and below 
this line. Because we assume spherical symmetry, we define 
the two caustics to be V„(V) = ±min{|t) up (r)|, |«down(r)|}- 
This choice limits the number of interlopers. The caustic 
amplitude is A K (r) — \V^(r) — V t r(r)]/2. Any realistic mod- 
els of galaxy clusters has dln^l/dlnr < 1/4. A value of 
this derivative much larger than 1/4 indicates that the algo- 
rithm has found an incorrect location of the caustics, which 
includes an excessive number of foreground and/or back- 
ground galaxies. Therefore, whenever dln^l/dlnr > ( at a 
given r, the algorithm replaces A(r) with a new value that 
yields dln^l/dlnr = 1/4. We choose £ = 2, rather than 
( — 1/4, to keep this constrain loose. 

We choose the parameter k, that determines the correct 
caustic location, as the root of the equation 



S(k) = (u csc ) k ,h -4{v ) = , 



(22) 



where {v eBC ) K ,i 



J Q R A 2 K (r)(p(r)dr / J <p(r)dr is the 
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parameter 


symbol 


default value 


Binary tree construction and candidate members 


galaxy mass 


in 


10 12 ft- 1 M s 


cluster threshold node 






substructure threshold node 


X2 




Caustic location 




rescaling 




25 


smoothing length 


h c 




threshold 


K 




derivative limit of A 


c 


2 



Mass profile 

filling factor Tp 0.7 

Table 1. List of the tunable parameters. 



caustic amplitude within the main group size R, ip(r) = 
J /q( r >' u )d«, and {v 2 } 1 ^ 2 is the velocity dispersion of the 
candidate cluster members with respect to the median red- 
shift. When S(k) — 0, the escape velocity inferred from the 
caustic amplitude equals the escape velocity of a system in 
dynamical equilibrium with a Maxwellian velocity distribu- 
tion within R. However, it is important to emphasize that 
the prescription for determining k works equally well for 
any cluster independently of its dynamical status within R. 
Therefore this prescription appears to be just a convenient 
recipe to locate the caustics properly and does not limit the 
caustic technique to systems in equilibrium within R. 

The outcome of the process described in the preceding 
paragraphs is outlined in Figure El where V±(r) are shown 
together with the expected caustic amplitude v CBCt i OB (r) — 
±y/2 | <p(r) | Jg{0) = ±y/(j) (r). Figure U shows a good 
agreement between the expected and estimated caustics; in 
Sect. 15.21 we will see that this is a general result. 



4.4 Tunable parameters 

The caustic technique requires two sets of parameters: a set 
for building the binary tree and identifying the cluster can- 
didate members and a set for the location of the caustics 
in the redshift diagram. Most parameters are automatically 
determined by the method with the prescriptions described 
above. A few remaining parameters, listed in Table [T] are 
chosen in input. In most cases, keeping these parameters to 
the default value returns accurate results. In fact, in Section 
16.21 we show how, in general, the improvement provided by 
a fine tuning of these parameters is modest. 



5 THE CAUSTIC TECHNIQUE AT WORK 

We now apply the caustic method to the mock redshift cata- 
logues described in Sect. 13.21 We describe how the technique 
recovers the centre and velocity dispersion of a cluster, and 
its escape velocity and mass profiles. Finally, we provide a 
simple recipe to estimate the uncertainties on these profiles. 



5.1 Global properties of the cluster 

In order to remove foreground and background large-scale 
structure from the redshift diagram, we remove all the par- 




0.0 



2 

r/r 



200 

Figure 10. Profiles of the ratio between the caustic ampli- 
tude A(r) and the l.o.s. component of the true escape velocity 
('"esc los( r )) = — 2 < / , ( r )/9(/3) = 4'i3( r )- The numerical gravitational 
potential profile <j>(r) is derived from the true mass distribution 
within r max = lO/i" 1 Mpc from the cluster centre: 50, 68, and 
90 percent of the profiles are within the upper and lower solid, 
dashed, and dotted curves, respectively. The solid squares show 
the median profile. The darkness of the shaded areas is propor- 
tional to the profile number density on the vertical axis. 



tides with l.o.s. velocity larger than 3.5 times the velocity 
dispersion of the candidate members. Nevertheless, when 
the cluster is embedded in a particularly crowded area of 
the sky, it can also happen that the main branch of the tree 
identifies a cluster different from the target cluster. These 
cases can happen, for example, when the target cluster is less 
rich than nearby systems; they are easily solved by limiting 
the input particle catalogue to a small enough region cen- 
tred on the target cluster. In our sample, the main branch of 
the binary tree identifies a cluster different from our target 
326 times out of 3000, in other words only 11 percent of the 
times. For the sake of simplicity, we limit our sample to the 
2674 clusters which are identified without further limitation 
of the cluster field of view. 

The input centre [a,8,cz] = [90° , 0°, 32000 kms -1 ], 
is recovered very well: the percentile ranges are a = 
[89.96,90.00,90.04]°, 6 = [-0.040,0.000,0.036]°, and cz = 
[31740,31998,32222] km s -1 . 

Figure[5]shows the agreement between the estimated ve- 
locity dispersion er C aus (the velocity dispersion of the galaxies 
hanging from the node xi) and the true velocity dispersion 
Otruej defined as the l.o.s. velocity dispersion of the parti- 
cles within a sphere of radius 3r2oo in real space: in 50 (95) 
percent of the systems the estimated velocity dispersion is 
within 5 (30) percent of the real one. 



5.2 The escape velocity profile 

The caustic amplitude is related to <j>p(r) (equation [2}, the 
mean component along the line of sight of the escape ve- 
locity, which is a combination of the gravitational poten- 
tial and the velocity anisotropy parameter. Figure [101 shows 
how well the caustic amplitude recovers </>g(r). On average, 
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Figure 9. Comparison between the true velocity dispersion ctrue and the velocity dispersion of the cluster candidate members <r C aus- 
The red solid line in the histogram in the right panel represents the median and the dash-dotted orange (dashed violet) lines represent 
the 50 (90) percent range. 



the potential is slightly underestimated in the central region 
(r < 0.2r2oo), but in the outer regions the agreement with 
the true escape velocity is remarkable. The slight systematic 
underestimate at the centre is due to the small number of 
particles at very small radii in the redshift diagram, that 
leads to an underestimate of the caustic amplitude. In the 
cluster outskirts, the location of the caustics is less accurate 
because the number of particles sampling the density distri- 
bution f q (r,v) is smaller. This fact increases the scatter of 
the profiles at large radii. We will see below that increasing 
the number of particles in the mock catalogues indeed re- 
duces the scatter (Figure [20)l . We remind, however, that the 
statistical significance of the scatter decreases at increasing 
radii: in fact, the number of clusters contributing to the me- 
dian profile decreases at large r, because we exclude those 
clusters that have a null caustic amplitude A(r) = beyond 
any given radius r. 

5.3 The mass profile 

DG97 show that the caustic amplitude can be related to the 
cumulative mass profile of the cluster by the relation 



GM(< r) 



T p {r)A 2 {r)dr 



The bottom right panel of Figure|4]shows the function Tp (r) 
in our simulations. At radii in the range ~ (0.5 — 4)r2oo, the 
average Tp{f) has a mild variation, between 0.5 and 0.8. 
This result led DG97 and D99 to assume Tp(r) = const 
tout court and assume that the mass profiles of real clusters 
can be estimated with the expression ©: 



GM(< r) 



Tp / A\r)dr 
Jo 



We can choose the correct value of the factor by con- 
sidering the contribution of the filling function J~p(r) in 
the integral of equation @. Figure [TT] shows (J~s(r)) = 
Jq Tp{x)Ax/r, where J-p(x) is the profile of each individual 
cluster. At radii larger than ~ 0.5r2oo, {J~p(r)) is basically 
constant and supports the validity of equation ©. 

We see that the most appropriate value is Tp = 0.7. 
This choice disagrees with the value Tp = 0.5 adopted by 




Figure 11. Profiles of the integral J Q Tp (x)dx/r described in the 
text; 90, 68, and 50 percent of the profiles are within the upper 
and lower dashed, solid and dotted curves. The central solid curve 
is the median profile. 



DG97 and D99. In this early work, the algorithm for the 
determination of the a plateau was less accurate than our 
algorithm here and systematically provided slightly larger 
caustic amplitudes. This overestimate was compensated by 
a smaller Tp, that, on turn, returned the correct mass pro- 
file, on average. Here, our improved algorithm appears to be 
more appropriate because it returns the correct <j>p(r) pro- 
file (Figure I10|l and, in order to estimate the correct mass 
profile, requires a value of Tp in agreement with what can 
be expected by inspecting Figure [TT] 

Figure [T2] shows that, on average, the mass profile is 
estimated at better than 10 percent at radii larger than ~ 
0.6r2oo- Clearly, at smaller radii, the assumption Tp(r) = 
const breaks down and the mass is severely overestimated. 

As already suggested by DG97, if we assume that the 
cluster is in virial equilibrium in the central region, we can 
use the virial theorem to estimate the mass there and limit 
the use of the caustic method to the cluster outskirts alone, 
where the equilibrium assumption does not hold. Here we 
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use the virial t heorem and the med ian and average mass 
estimators from iHeisler et all (|l985h to estimate the mass 
within aR, where R is the mean clustrocentric separation 
of the candidate cluster members from the binary tree (see 
Section 14-3 P and a is a free parameter. In our sample, the 
percentile range is R = [0.50, 1.23, 1.68]/i _1 Mpc. We com- 
pute the ratio between the estimated mass and the true mass 
for different values of a. We find that the best estimates are 
obtained when a = 0.7. In this case, the ratio between the 
estimated and true masses is, on average, 1.03 for the virial 
theorem, 1.30 and 1.49 for the median and average mass 
estimators, respectively. Different values of a yield worse 
mass estimates. This result indicates that the radius 0.7R 
generally contains the cluster region in approximate virial 
equilibrium. 

When we estimate the mass with the virial theorem 
within 0.7R and use equation §5} with 0.7R as the lower 
limit of the integral, we still obtain a very good estimate of 
the real mass (Figure [T3J. 



5.4 The gravitational potential profile 

Within r — l/i -1 Mpc, we fit the caustic mass profile of each 
individual cluster with an NFW profile. The scale factor 
?"200 derived from this fit agrees with the true 7-200 derived 
in 3D; the percentile range of their ratio is [0.84, 1.03, 1.33]. 
Consequently, for M200, which is proportional to rfoo, we 
have [0.59,1.10,2.36] (Figure OH) . 

However, the concentration parameter c derived from 
the NFW fit to the caustic mass profile tends to overesti- 
mate by 36 percent, on average, the concentration param- 
eter of the NFW fit to the 3D mass profile: in fact, the 
caustic c has percentile range [3.59, 5.92, 10.23], whereas we 
find [2.06,4.62,6.50] in 3D; the ratio of these two parame- 
ters has percentile range [0.73, 1.36,3.26]. This result is an 
obvious consequence of the mass overestimate at small radii 
(Figure EJ). 

The NFW fit parameters can be used to derive the grav- 



itational potential profile of the cluster. The left panel of 
Figure [15] shows that the estimated gravitational potential 
profile returns the true profile within 10 percent, on aver- 
age, and, in 50 percent of the cases, the estimated profile 
is within 25 percent from the real one out to 4r2oo- This 
result is impressive and derives from the mild radial varia- 
tion of the functions J-(r) and g(/3) in hierarchical clustering 
models. 

The right panel of Figure [TS] shows the ratio between 
the caustic mass profile derived from the NFW fit and the 
true mass profile. The agreement is again good, within 20 
percent, although not as good as in Figure [T2l because the 
caustic estimate is biased low due to the overestimate of c. 

5.5 Uncertainty estimate 

Figures[lO]and[T2]show the spread derived from the distribu- 
tion of the profiles of the individual clusters. This spread is 
mostly due to projection effects, as we will see below (Sect. 
16 . 3 P - We remind here the simple recipe suggested in D99 to 
estimate this spread when observing a real cluster of galax- 
ies. 

The uncertainty in the measured value of A(r) depends 
on the number of galaxies contributing to the determination 
of A(r); we thus define the relative error 

SA(r)/A(r) = k/ max{f q (r,v)}, (23) 

where the maximum value of f q (r,v) is along the v-axis at 
fixed r. The resulting error on the cumulative mass profile 
is 

5M i = ^ ^mjSAir^/A^)], (24) 

where raj is the mass of the shell [rj-i, rj] and i is the index 
of the most external shell. 

The shaded bands in Figure [16] are the median spreads 
derived from the above equations. The median uncertainty 
on the escape velocity profile increases up to r ~ 1.5r2oo 
and remains within ~ 20 — 30 percent at larger radii. The 
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Figure 14. Distribution of the ratio between the parameters of the best NFW fit to the caustic mass profile and the true M200 and T200 
(left and central panels) and distribution of the concentration parameter of the best NFW fit to the caustic mass profile (right panel). 
The dashed lines indicate the median while the solid lines show the 90 percent range. 
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Figure 15. Left panel: Profiles of the ratio between the NFW gravitational potential profile, with parameters derived from the NFW 
best fit to the caustic mass profile, and the numerical gravitational potential profile derived from the true mass distribution within 
r maI = 10h~ 1 Mpc from the cluster centre. Right panel: Profiles of the ratio between the NFW best fit to the caustic mass profile, 
and the true mass profile. In both panels, 90 and 50 percent of the profiles are within the upper and lower dotted and solid curves, 
respectively. The central solid curves are the median profiles. For each cluster, the NFW fit to the mass profile is only performed to the 
mass distribution within lh~ 1 Mpc. 
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Figure 16. The shaded areas show the median uncertainties estimated from equations H23H and (124[ i. For comparison, we also show the 
90 (dotted lines), 68 (dashed lines), and 50 (solid lines) percent ranges of the profiles from Figure llOl The solid squares show the median 
profiles. 
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Figure 17. Profiles of the ratio between the caustic and the true 
mass profile when the centre is forced to be the true one. The 
lines and shaded areas are as in Figure fl6l 

mass profile has a slightly larger uncertainties but remains 
below ~ 30 — 40 percent out to r ~ 3r2oo- At radii smaller 
than ~ 0.6r2oo the assumption of a constant filling function 
J-p(r) breaks down, as we mentioned above, but this mass 
overestimate is systematic and can be easily taken into ac- 
count. 

Figure [16] shows that our recipes for the estimate of 
the uncertainties yield values in good agreement with the 
50 percent range of the profiles. Therefore, when applied to 
real clusters, this algorithm provides uncertainties with 50 
percent confidence levels. 



6 SYSTEMATICS 

In this section we investigate the systematic errors affecting 
the estimate of the escape velocity and mass profiles: the 
systematic errors can originate from the centre location, the 
parameters of the algorithm, and the projection effects. We 
also show how the accuracy of the estimate depends on the 
number of objects in the redshift catalogue. 

6.1 Centre identification 

In Sect. 15. ll we show how the cluster centre generally is well 
recovered. The deviations are smaller than 0.07° on the sky 
and 250 km s _1 along the line of sight in 90 percent of the 
clusters. These deviations produce negligible effects on the 
mass profile estimate. Figure [T7] shows the ratio between 
the caustic mass profile and the true mass profile, when the 
centre of the cluster is imposed to be the true one. A com- 
parison with Figure [T2l shows that the average profile is still 
correct at radii larger than 0.6r2oo: forcing the method to 
use the correct centre only slightly improves the estimate 
and reduces the scatter at radii larger than ~ 2.5r2oo- It is 
worth pointing out that forcing the code to use the correct 
centre alters the quantities involved in the construction of 
the redshift diagram, specifically we recompute the veloc- 
ity dispersion and the cluster size R of the candidate cluster 
members with respect to the new centre. We then derive the 



new redshift diagram and locate the caustics. Despite these 
variations, the mass profile estimate does not substantially 
change. This result confirms the robustness of the method. 



6.2 Tuning the parameters 

The algorithm allows the user to choose some parame- 
ters, namely the thresholds of the binary tree, the optimal 
smoothing length h c for the galaxy density distribution in 
the redshift diagram, and the threshold k. This freedom is 
necessary when the target cluster is in particularly crowded 
regions, or the galaxy sampling is too sparse. In these cases, 
the algorithm is either unable to locate the caustics or it 
returns unrealistic caustic amplitudes. If one is interested 
in the average profiles of a large cluster sample, tuning the 
caustic parameters for each individual cluster can be very 
time consuming. To quantify the impact of this possible fine 
tuning, for each of our 100 simulated clusters, we randomly 
draw one of the 30 projections and tune the input parame- 
ters by hand until the caustic location appears to be close 
enough to where we might expect them to be by eye. It turns 
out that most clusters need a rather minor fine tuning, as 
demonstrated by the final result shown in Figure 1181 The 
left-hand panel shows the caustic amplitude when we ap- 
ply a fine tuning to the parameters; the right-hand panel is 
for the same sample without fine tuning. Clearly, the me- 
dian profile remains unchanged, but the scatter is slightly 
reduced. 



6.3 Projection effects 

The analysis provided above clearly shows that the uncer- 
tainties on the cluster centre determination and the freedom 
on the algorithm parameters are not responsible for most of 
the spread of the escape velocity and mass profiles. 

This spread originates from the assumption of spher- 
ical symmetry. In hierarchical clustering, this assumption 
does not hold in general, but observationally our informa- 
tion is limited to the galaxy distribution on the sky alone, 
although various techniques can in principle provide infor- 
mation on the 3D shape of the cluster (e.g. IZaroubi et al.l 
l200ll : lAmegTio et al ] l2009h . 

To show the impact of the projection effects on the caus- 
tic method, we take our 100 clusters and plot the escape 
velocity profiles derived from each of the 30 lines of sight. 
Figure [19] shows four randomly chosen clusters as examples. 
The caustic technique returns a median profile systemati- 
cally large for the cluster in the top-right panel. In the other 
three cases, however, the median profile is within 30 percent 
of the correct one out to ~ 3r2oo • 

The relevant result of this test is the fact that the spread 
due to the different lines of sight is comparable to the spread 
of the entire sample (shaded area) taken from Figure 1101 
This result clearly indicates that the projection effects are 
the major responsible for the systematic uncertainties of the 
caustic technique and further refinements of the technique, 
which still assume spherical symmetry, appear to be unable 
to improve the mass estimate. 
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Figure 18. Profiles of the ratio between the caustic and the true escape velocity profiles when the parameters are finely tuned (left 
panel) and from the automatic procedure (right panel) in a subsample of 100 mock catalogues: 50, 68, and 90 percent of the profiles are 
within the upper and lower solid, dashed, and dotted curves, respectively. The solid squares show the median profiles. 




Figure 19. Profiles of the ratio between the caustic amplitude A(r) and the l.o.s. component of the true escape velocity {v^ i os (?*)) = 
— 2</>(r) / g(/3) = 4>{i{ r ) f° r f° ur randomly chosen clusters. The numerical gravitational potential profile (j>(r) is derived from the true mass 
distribution within r max = 10h~ 1 Mpc from the cluster centre: 50, 68, and 90 percent of the profiles derived from 30 different l.o.s. are 
within the upper and lower solid, dashed, and dotted curves, respectively. The solid squares show the median profiles. The shaded areas 
are taken from Figure ITOl 
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Figure 20. Ratio between the estimated and true escape velocity profiles for subsamples with a given number of particles within 3r"200 as 
indicated in each panel; 50, 68, and 90 percent of the profiles are within the upper and lower solid, dashed, and dotted curves, respectively. 
The solid squares show the median profiles. 



6.4 Dependence of the profiles on the number of 
particles 

To return a reliable estimate of the escape velocity and mass 
profiles, the caustic technique requires a sufficient number 
of redshifts sampling the cluster velocity field. 

To quantify this number, we stack our 3000 mock cata- 
logues and apply the caustic technique to the stacked clus- 
ter with an increasing number of particles in the catalogue. 
With this procedure we average over the different lines of 



sight and the final stacked cluster is very close to be spheri- 
cally symmetric. We randomly choose particles in the cata- 
logue until we obtain a given number N of particles within 
3r2oo from the cluster centre in real space. We consider cata- 
logues with TV = [100, 200, 350, 500, 1000, 2000]. We compile 
100 different catalogues for each value of TV. The catalogues 
with the same TV clearly have different total numbers of par- 
ticles in the catalogue. Table [2] lists the percentile ranges of 
the total number of particles in the catalogues, which, as 
usual, have a field of view 12 x 12 h~ 2 Mpc 2 at the clus- 
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ter redshift. Table [2] also lists the total number of particles 
in a field of view of 30' x 30' (2.46 x 2.46 h~ 2 Mpc 2 ) and 
|Aui os | ^ 2000 km s" 1 centred on the cluster. 

Figure I2TJ1 shows the escape velocity profiles for the six 
different sets of catalogues. The spread decreases with the in- 
creasing number of particles, as expected. However, the pro- 
file tends to be overestimated at large radii and with richer 
catalogues. This result origi nates from the fact that, as i t 
was already noticed by D99 , ISchmalzing fc Diaferid |2000), 
and ICasagrande fc Diaferid (|2006ft . mock redshift surveys 
extracted from Af-body simulations show large-scale struc- 
tures which are less sharp than in the real universe. There- 
fore, when we increase the number of particles in our mock 
catalogues, we start sampling the surrounding structure of 
the cluster in real space, the expected trumpet shape of the 
particle density distribution in the redshift diagram becomes 
fuzzier, and the caustic algorithm has increasing difficulty at 
locating the causti cs properly. This proble m is not present 
in real surveys (e.g. lRines fc D iafcrio 20061') . At these scales, 
the dark matter particles, that we use here, and the galax- 
ies might not share the same distribution in the redshift 
diagram and an appropriate galaxy formation algorithm in 
mock surveys might provide sharper large-scale structures. 
At any rate, Figure [2U1 shows that as few as ~ 200 galaxies 
within 3r2oo, i-e. a few tens of redshifts per squared comov- 
ing megaparsec within the cluster, are sufficient to have a 
reasonably accurate estimate of the escape velocity profile 
out to ~ 4r2oo- 



7 DISCUSSION AND CONCLUSION 

We use a large sample of galaxy clusters extracted from 
a cosmological TV-body simulation to show that, by using 
galaxy redshifts alone, the caustic technique measures (1) 
the escape velocity profiles of individual galaxy clusters with 
1-cr uncertainty between ~ 20 and ~ 50 percent when the 
clustrocentric distances increase from ~ 0.1r2oo to 4r2oo, and 
(2) the cluster mass profiles with ~ 50 percent 1-cr uncer- 
tainty at clustrocentric distances in the range ~ (0.6— 4)r2oo- 
The technique returns reliable estimates when a few tens of 
redshifts per squared comoving megaparsec within the clus- 
ter are available. For sparser samples, the technique can be 
applied to synthetic clusters obtained by stacking individ- 
ual clusters. The technique relies on the kinematic proper- 
ties of galaxy clusters in hierarchical clustering scenarios of 
structure formation and assumes spherical symmetry. This 
assumption alone is responsible for most of the estimated 
uncertainty. 

The caustic amplitude measures the profile of the es- 
cape velocity along the line of sight, which is a combina- 
tion of the gravitational potential profile and the velocity 
anisotropy parameter profile /3(r). There are extensive at- 
tempts to estimate f5(r) and mass profiles by modelling the 
density of galaxies in the redshift diagrams with sensible 
distri bution functions l|Woitak et al.ll2009l ; IWoitak fc Lokasl 
l2010l ). These methods assume dynamical equilibrium and 
the caustic measure can thus provide a robust consistency 
che ck of these estimates of f3(r) within r^oo- 

ICupani et alj (|2008T ) have proposed a new approach, 
based on the spherical collapse model, to estimate the mass 
profile of galaxy clusters beyond their virial radius. The 



spherical collapse model poorly describes the formation of 
individual clusters in hi e rarchi cal clustering scenarios. Nev- 
ertheless, ICupani et al. (2008) 's method seems to provide 
cluster masses with 50 percent uncertainty on average, once 
the mean density of the universe Qo is set and an appropriate 
subsample o f galaxies is chosen based on their phase-space 

coordinates l|Cupani et al.ll2010h . 

In addition to the caustic method and the lCupani et all 

(2008)'s technique, the only method available to estimate 
the mass of clusters beyond their virial radius is gravita- 
tional lensing. However, the lensing signal is strong enough 
only when the cluster is within the redshift range z ~ 0.1 — 1, 
whereas, in principle, the caustic technique can be applied 
to clusters at any redshift and it is only limited by the ob- 
serving time required to measure a large enough number of 
galaxy spectra. The caustic technique relies on the assump- 
tion that randomly chosen galaxies are fair tracers of the 
velocity field, whereas gravitational lensing provides a di- 
rect measure of the projected mass distribution. All meth- 
ods based on optical data relies on the absence of veloc- 
ity bias between gala xies and dark matter, and both N- 
body simulations (e. g. iDiaferio et al.l l200ll ; iGill et all |2004| ; 
Diemand et al. ||2004l ; lGill et al.ll2005l ) and observations (e.g. 



Rines et al.l 12008^ indicate that this assumption is indeed 



reasonable. Selected galaxy samples do however have ve- 
locity biases: in fact, elliptical galaxies, disk galaxies, and 
galaxies showing signs of interactions can have substan- 
tially diff erent velocity distributions bot h in observations 
(e.g. iBiviano fc Katgert]|2q04l ; [Mos3|2006| ') and in simulated 
clusters (e.g. |Piaferidll999T l. 

The caustic technique estimates the escape velocity 
from a cluster. Therefore, it can also be used to identify the 
galaxy members of the cluster. This argument is not limited 
to clusters, and the caustic technique was indee d applied to 
remo ving stellar interlopers in dwarf g alaxies dSerra et al.l 
l20ld ) and the Milky Way stellar halo (|Brown et alj|2010h . 
The technique arranges the galaxies in the cluster field in a 
binary tree and a by-product of this analysis is the substruc- 
ture identification. Here, we have focussed on the estimate 
of the escape velocity and mass profiles of galaxy clusters. 
We will investigate the use of the caustic technique to iden- 
tify members and substructure of self-gravitating systems in 
future work. 

With the advent of the future generation of large spec- 
troscopic surveys, both from ground (e.g., BigBOSjjf]) and 
from space (e.g., EUCLIE0), the application of the caustic 
technique can represent an extremely valuable tool to esti- 
mate mass profiles and galaxy membership over a large en- 
semble of galaxy clusters out to large radii from the cluster 
centre, thus providing important information on the cosmo- 
logical assembly of dark matter halos and on the galaxy- 
environment connection. 
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Table 2. Number of particles in the catalogues compiled from the stacked cluster. The FOV columns list the numbers of particles within 
the field of view 30' X 30' (2.46 X 2.46 h' 2 Mpc 2 ) and \Av ios \ sC 2000 km s" 1 centred on the cluster. 
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